1 /* Copyright 2002-2016 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.semianalytical.dsst.forces;
18
19 import java.io.NotSerializableException;
20 import java.io.Serializable;
21 import java.util.ArrayList;
22 import java.util.HashMap;
23 import java.util.List;
24 import java.util.Map;
25 import java.util.Set;
26 import java.util.SortedSet;
27
28 import org.hipparchus.analysis.UnivariateVectorFunction;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.hipparchus.util.FastMath;
31 import org.orekit.attitudes.Attitude;
32 import org.orekit.attitudes.AttitudeProvider;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitExceptionWrapper;
35 import org.orekit.forces.ForceModel;
36 import org.orekit.frames.Frame;
37 import org.orekit.orbits.EquinoctialOrbit;
38 import org.orekit.orbits.Orbit;
39 import org.orekit.orbits.PositionAngle;
40 import org.orekit.propagation.SpacecraftState;
41 import org.orekit.propagation.numerical.TimeDerivativesEquations;
42 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
43 import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
44 import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
45 import org.orekit.time.AbsoluteDate;
46 import org.orekit.utils.TimeSpanMap;
47
48 /** Common handling of {@link DSSTForceModel} methods for Gaussian contributions to DSST propagation.
49 * <p>
50 * This abstract class allows to provide easily a subset of {@link DSSTForceModel} methods
51 * for specific Gaussian contributions.
52 * </p><p>
53 * This class implements the notion of numerical averaging of the DSST theory.
54 * Numerical averaging is mainly used for non-conservative disturbing forces such as
55 * atmospheric drag and solar radiation pressure.
56 * </p><p>
57 * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = δa<sub>i</sub>/δv . q<br>
58 * where:
59 * <ul>
60 * <li>a<sub>i</sub> are the six equinoctial elements</li>
61 * <li>v is the velocity vector</li>
62 * <li>q is the perturbing acceleration due to the considered force</li>
63 * </ul>
64 *
65 * <p> The averaging process and other considerations lead to integrate this contribution
66 * over the true longitude L possibly taking into account some limits.
67 *
68 * <p> To create a numerically averaged contribution, one needs only to provide a
69 * {@link ForceModel} and to implement in the derived class the method:
70 * {@link #getLLimits(SpacecraftState)}.
71 * </p>
72 * @author Pascal Parraud
73 */
74 public abstract class AbstractGaussianContribution implements DSSTForceModel {
75
76 /** Available orders for Gauss quadrature. */
77 private static final int[] GAUSS_ORDER = {12, 16, 20, 24, 32, 40, 48};
78
79 /** Max rank in Gauss quadrature orders array. */
80 private static final int MAX_ORDER_RANK = GAUSS_ORDER.length - 1;
81
82 /** Number of points for interpolation. */
83 private static final int INTERPOLATION_POINTS = 3;
84
85 /** Maximum value for j index. */
86 private static final int JMAX = 12;
87
88 /** Retrograde factor I.
89 * <p>
90 * DSST model needs equinoctial orbit as internal representation.
91 * Classical equinoctial elements have discontinuities when inclination
92 * is close to zero. In this representation, I = +1. <br>
93 * To avoid this discontinuity, another representation exists and equinoctial
94 * elements can be expressed in a different way, called "retrograde" orbit.
95 * This implies I = -1. <br>
96 * As Orekit doesn't implement the retrograde orbit, I is always set to +1.
97 * But for the sake of consistency with the theory, the retrograde factor
98 * has been kept in the formulas.
99 * </p>
100 */
101 private static final int I = 1;
102
103 // CHECKSTYLE: stop VisibilityModifierCheck
104
105 /** a. */
106 protected double a;
107 /** e<sub>x</sub>. */
108 protected double k;
109 /** e<sub>y</sub>. */
110 protected double h;
111 /** h<sub>x</sub>. */
112 protected double q;
113 /** h<sub>y</sub>. */
114 protected double p;
115
116 /** Eccentricity. */
117 protected double ecc;
118
119 /** Kepler mean motion: n = sqrt(μ / a³). */
120 protected double n;
121
122 /** Mean longitude. */
123 protected double lm;
124
125 /** Equinoctial frame f vector. */
126 protected Vector3D f;
127 /** Equinoctial frame g vector. */
128 protected Vector3D g;
129 /** Equinoctial frame w vector. */
130 protected Vector3D w;
131
132 /** A = sqrt(μ * a). */
133 protected double A;
134 /** B = sqrt(1 - h² - k²). */
135 protected double B;
136 /** C = 1 + p² + q². */
137 protected double C;
138
139 /** 2 / (n² * a) . */
140 protected double ton2a;
141 /** 1 / A .*/
142 protected double ooA;
143 /** 1 / (A * B) .*/
144 protected double ooAB;
145 /** C / (2 * A * B) .*/
146 protected double co2AB;
147 /** 1 / (1 + B) .*/
148 protected double ooBpo;
149 /** 1 / μ .*/
150 protected double ooMu;
151 /** μ .*/
152 protected double mu;
153
154 // CHECKSTYLE: resume VisibilityModifierCheck
155
156 /** Contribution to be numerically averaged. */
157 private final ForceModel contribution;
158
159 /** Gauss integrator. */
160 private final double threshold;
161
162 /** Gauss integrator. */
163 private GaussQuadrature integrator;
164
165 /** Flag for Gauss order computation. */
166 private boolean isDirty;
167
168 /** Attitude provider. */
169 private AttitudeProvider attitudeProvider;
170
171 /** Prefix for coefficients keys. */
172 private final String coefficientsKeyPrefix;
173
174 /** Short period terms. */
175 private GaussianShortPeriodicCoefficients gaussianSPCoefs;
176
177 /** Build a new instance.
178 * @param coefficientsKeyPrefix prefix for coefficients keys
179 * @param threshold tolerance for the choice of the Gauss quadrature order
180 * @param contribution the {@link ForceModel} to be numerically averaged
181 */
182 protected AbstractGaussianContribution(final String coefficientsKeyPrefix,
183 final double threshold,
184 final ForceModel contribution) {
185 this.coefficientsKeyPrefix = coefficientsKeyPrefix;
186 this.contribution = contribution;
187 this.threshold = threshold;
188 this.integrator = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
189 this.isDirty = true;
190 }
191
192 /** {@inheritDoc} */
193 @Override
194 public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly) {
195
196 final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
197 gaussianSPCoefs = new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix,
198 JMAX, INTERPOLATION_POINTS,
199 new TimeSpanMap<Slot>(new Slot(JMAX, INTERPOLATION_POINTS)));
200 list.add(gaussianSPCoefs);
201 return list;
202
203 }
204
205 /** {@inheritDoc} */
206 @Override
207 public void initializeStep(final AuxiliaryElements aux)
208 throws OrekitException {
209
210 // Equinoctial elements
211 a = aux.getSma();
212 k = aux.getK();
213 h = aux.getH();
214 q = aux.getQ();
215 p = aux.getP();
216
217 // Eccentricity
218 ecc = aux.getEcc();
219
220 // Equinoctial coefficients
221 A = aux.getA();
222 B = aux.getB();
223 C = aux.getC();
224
225 // Equinoctial frame vectors
226 f = aux.getVectorF();
227 g = aux.getVectorG();
228 w = aux.getVectorW();
229
230 // Kepler mean motion
231 n = aux.getMeanMotion();
232
233 // Mean longitude
234 lm = aux.getLM();
235
236 // 1 / A
237 ooA = 1. / A;
238 // 1 / AB
239 ooAB = ooA / B;
240 // C / 2AB
241 co2AB = C * ooAB / 2.;
242 // 1 / (1 + B)
243 ooBpo = 1. / (1. + B);
244 // 2 / (n² * a)
245 ton2a = 2. / (n * n * a);
246 // mu
247 mu = aux.getMu();
248 // 1 / mu
249 ooMu = 1. / mu;
250 }
251
252 /** {@inheritDoc} */
253 @Override
254 public double[] getMeanElementRate(final SpacecraftState state) throws OrekitException {
255
256 double[] meanElementRate = new double[6];
257 // Computes the limits for the integral
258 final double[] ll = getLLimits(state);
259 // Computes integrated mean element rates if Llow < Lhigh
260 if (ll[0] < ll[1]) {
261 meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1]);
262 if (isDirty) {
263 boolean next = true;
264 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
265 final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1]);
266 if (getRatesDiff(meanElementRate, meanRates) < threshold) {
267 integrator = new GaussQuadrature(GAUSS_ORDER[i]);
268 next = false;
269 }
270 }
271 isDirty = false;
272 }
273 }
274 return meanElementRate;
275 }
276
277 /** Compute the acceleration due to the non conservative perturbing force.
278 *
279 * @param state current state information: date, kinematics, attitude
280 * @return the perturbing acceleration
281 * @exception OrekitException if some specific error occurs
282 */
283 protected Vector3D getAcceleration(final SpacecraftState state)
284 throws OrekitException {
285 final AccelerationRetriever retriever = new AccelerationRetriever(state);
286 contribution.addContribution(state, retriever);
287
288 return retriever.getAcceleration();
289 }
290
291 /** Compute the limits in L, the true longitude, for integration.
292 *
293 * @param state current state information: date, kinematics, attitude
294 * @return the integration limits in L
295 * @exception OrekitException if some specific error occurs
296 */
297 protected abstract double[] getLLimits(final SpacecraftState state) throws OrekitException;
298
299 /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
300 *
301 * @param state current state
302 * @param gauss Gauss quadrature
303 * @param low lower bound of the integral interval
304 * @param high upper bound of the integral interval
305 * @return the mean element rates
306 * @throws OrekitException if some specific error occurs
307 */
308 private double[] getMeanElementRate(final SpacecraftState state,
309 final GaussQuadrature gauss,
310 final double low,
311 final double high) throws OrekitException {
312 final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0), low, high);
313 // Constant multiplier for integral
314 final double coef = 1. / (2. * FastMath.PI * B);
315 // Corrects mean element rates
316 for (int i = 0; i < 6; i++) {
317 meanElementRate[i] *= coef;
318 }
319 return meanElementRate;
320 }
321
322 /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
323 *
324 * @param meanRef reference rates
325 * @param meanCur current rates
326 * @return estimated magnitude of weighted differences
327 */
328 private double getRatesDiff(final double[] meanRef, final double[] meanCur) {
329 double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / a;
330 // Corrects mean element rates
331 for (int i = 1; i < meanRef.length; i++) {
332 final double diff = FastMath.abs(meanRef[i] - meanCur[i]);
333 if (maxDiff < diff) maxDiff = diff;
334 }
335 return maxDiff;
336 }
337
338 /** {@inheritDoc} */
339 @Override
340 public void registerAttitudeProvider(final AttitudeProvider provider) {
341 this.attitudeProvider = provider;
342 }
343
344 /** {@inheritDoc} */
345 @Override
346 public void updateShortPeriodTerms(final SpacecraftState ... meanStates)
347 throws OrekitException {
348
349 final Slot slot = gaussianSPCoefs.createSlot(meanStates);
350 for (final SpacecraftState meanState : meanStates) {
351 initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));
352 final double[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate());
353 final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(meanState, JMAX);
354 final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, JMAX);
355 gaussianSPCoefs.computeCoefficients(meanState, slot, fourierCjSj, uijvij, n, a);
356 }
357
358 }
359
360 /**
361 * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
362 * <p>
363 * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
364 * ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
365 * σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
366 * </p>
367 * @param date current date
368 * @return computed coefficients
369 */
370 private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date) {
371 final double[][] currentRhoSigmaj = new double[2][3 * JMAX + 1];
372 final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
373 final double b = 1. / (1 + B);
374
375 // (-b)<sup>j</sup>
376 double mbtj = 1;
377
378 for (int j = 1; j <= 3 * JMAX; j++) {
379
380 //Compute current rho and sigma;
381 mbtj *= -b;
382 final double coef = (1 + j * B) * mbtj;
383 currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
384 currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
385 }
386 return currentRhoSigmaj;
387 }
388
389 /** Internal class for retrieving acceleration from a {@link ForceModel}. */
390 private static class AccelerationRetriever implements TimeDerivativesEquations {
391
392 /** acceleration vector. */
393 private Vector3D acceleration;
394
395 /** state. */
396 private final SpacecraftState state;
397
398 /** Simple constructor.
399 * @param state input state
400 */
401 AccelerationRetriever(final SpacecraftState state) {
402 this.acceleration = Vector3D.ZERO;
403 this.state = state;
404 }
405
406 /** {@inheritDoc} */
407 @Override
408 public void addKeplerContribution(final double mu) {
409 }
410
411 /** {@inheritDoc} */
412 @Override
413 public void addXYZAcceleration(final double x, final double y, final double z) {
414 acceleration = new Vector3D(x, y, z);
415 }
416
417 /** {@inheritDoc} */
418 @Override
419 public void addAcceleration(final Vector3D gamma, final Frame frame)
420 throws OrekitException {
421 acceleration = frame.getTransformTo(state.getFrame(),
422 state.getDate()).transformVector(gamma);
423 }
424
425 /** {@inheritDoc} */
426 @Override
427 public void addMassDerivative(final double q) {
428 }
429
430 /** Get the acceleration vector.
431 * @return acceleration vector
432 */
433 public Vector3D getAcceleration() {
434 return acceleration;
435 }
436
437 }
438
439 /** Internal class for numerical quadrature. */
440 private class IntegrableFunction implements UnivariateVectorFunction {
441
442 /** Current state. */
443 private final SpacecraftState state;
444
445 /** Signal that this class is used to compute the values required by the mean element variations
446 * or by the short periodic element variations. */
447 private final boolean meanMode;
448
449 /** The j index.
450 * <p>
451 * Used only for short periodic variation. Ignored for mean elements variation.
452 * </p> */
453 private final int j;
454
455 /** Build a new instance.
456 * @param state current state information: date, kinematics, attitude
457 * @param meanMode if true return the value associated to the mean elements variation,
458 * if false return the values associated to the short periodic elements variation
459 * @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
460 */
461 IntegrableFunction(final SpacecraftState state, final boolean meanMode, final int j) {
462 this.state = state;
463 this.meanMode = meanMode;
464 this.j = j;
465 }
466
467 /** {@inheritDoc} */
468 @Override
469 public double[] value(final double x) {
470
471 //Compute the time difference from the true longitude difference
472 final double shiftedLm = trueToMean(x);
473 final double dLm = shiftedLm - lm;
474 final double dt = dLm / n;
475
476 final double cosL = FastMath.cos(x);
477 final double sinL = FastMath.sin(x);
478 final double roa = B * B / (1. + h * sinL + k * cosL);
479 final double roa2 = roa * roa;
480 final double r = a * roa;
481 final double X = r * cosL;
482 final double Y = r * sinL;
483 final double naob = n * a / B;
484 final double Xdot = -naob * (h + sinL);
485 final double Ydot = naob * (k + cosL);
486 final Vector3D vel = new Vector3D(Xdot, f, Ydot, g);
487
488 // Compute acceleration
489 Vector3D acc = Vector3D.ZERO;
490 try {
491
492 // shift the orbit to dt
493 final Orbit shiftedOrbit = state.getOrbit().shiftedBy(dt);
494
495 // Recompose an orbit with time held fixed to be compliant with DSST theory
496 final Orbit recomposedOrbit =
497 new EquinoctialOrbit(shiftedOrbit.getA(),
498 shiftedOrbit.getEquinoctialEx(),
499 shiftedOrbit.getEquinoctialEy(),
500 shiftedOrbit.getHx(),
501 shiftedOrbit.getHy(),
502 shiftedOrbit.getLv(),
503 PositionAngle.TRUE,
504 shiftedOrbit.getFrame(),
505 state.getDate(),
506 shiftedOrbit.getMu());
507
508 // Get the corresponding attitude
509 final Attitude recomposedAttitude =
510 attitudeProvider.getAttitude(recomposedOrbit,
511 recomposedOrbit.getDate(),
512 recomposedOrbit.getFrame());
513
514 // create shifted SpacecraftState with attitude at specified time
515 final SpacecraftState shiftedState =
516 new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());
517
518 acc = getAcceleration(shiftedState);
519
520 } catch (OrekitException oe) {
521 throw new OrekitExceptionWrapper(oe);
522 }
523 //Compute the derivatives of the elements by the speed
524 final double[] deriv = new double[6];
525 // da/dv
526 deriv[0] = getAoV(vel).dotProduct(acc);
527 // dex/dv
528 deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
529 // dey/dv
530 deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
531 // dhx/dv
532 deriv[3] = getQoV(X).dotProduct(acc);
533 // dhy/dv
534 deriv[4] = getPoV(Y).dotProduct(acc);
535 // dλ/dv
536 deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);
537
538 // Compute mean elements rates
539 double[] val = null;
540 if (meanMode) {
541 val = new double[6];
542 for (int i = 0; i < 6; i++) {
543 // da<sub>i</sub>/dt
544 val[i] = roa2 * deriv[i];
545 }
546 } else {
547 val = new double[12];
548 //Compute cos(j*L) and sin(j*L);
549 final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
550 final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);
551
552 for (int i = 0; i < 6; i++) {
553 // da<sub>i</sub>/dv * cos(jL)
554 val[i] = cosjL * deriv[i];
555 // da<sub>i</sub>/dv * sin(jL)
556 val[i + 6] = sinjL * deriv[i];
557 }
558 }
559 return val;
560 }
561
562 /** Converts true longitude to eccentric longitude.
563 * @param lv True longitude
564 * @return Eccentric longitude
565 */
566 private double trueToEccentric (final double lv) {
567 final double cosLv = FastMath.cos(lv);
568 final double sinLv = FastMath.sin(lv);
569 final double num = h * cosLv - k * sinLv;
570 final double den = B + 1 + k * cosLv + h * sinLv;
571 return lv + 2 * FastMath.atan(num / den);
572 }
573
574 /** Converts eccentric longitude to mean longitude.
575 * @param le Eccentric longitude
576 * @return Mean longitude
577 */
578 private double eccentricToMean (final double le) {
579 return le - k * FastMath.sin(le) + h * FastMath.cos(le);
580 }
581
582 /** Converts true longitude to mean longitude.
583 * @param lv True longitude
584 * @return Eccentric longitude
585 */
586 private double trueToMean (final double lv) {
587 return eccentricToMean(trueToEccentric(lv));
588 }
589
590 /** Compute δa/δv.
591 * @param vel satellite velocity
592 * @return δa/δv
593 */
594 private Vector3D getAoV(final Vector3D vel) {
595 return new Vector3D(ton2a, vel);
596 }
597
598 /** Compute δh/δv.
599 * @param X satellite position component along f, equinoctial reference frame 1st vector
600 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
601 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
602 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
603 * @return δh/δv
604 */
605 private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
606 final double kf = (2. * Xdot * Y - X * Ydot) * ooMu;
607 final double kg = X * Xdot * ooMu;
608 final double kw = k * (I * q * Y - p * X) * ooAB;
609 return new Vector3D(kf, f, -kg, g, kw, w);
610 }
611
612 /** Compute δk/δv.
613 * @param X satellite position component along f, equinoctial reference frame 1st vector
614 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
615 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
616 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
617 * @return δk/δv
618 */
619 private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
620 final double kf = Y * Ydot * ooMu;
621 final double kg = (2. * X * Ydot - Xdot * Y) * ooMu;
622 final double kw = h * (I * q * Y - p * X) * ooAB;
623 return new Vector3D(-kf, f, kg, g, -kw, w);
624 }
625
626 /** Compute δp/δv.
627 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
628 * @return δp/δv
629 */
630 private Vector3D getPoV(final double Y) {
631 return new Vector3D(co2AB * Y, w);
632 }
633
634 /** Compute δq/δv.
635 * @param X satellite position component along f, equinoctial reference frame 1st vector
636 * @return δq/δv
637 */
638 private Vector3D getQoV(final double X) {
639 return new Vector3D(I * co2AB * X, w);
640 }
641
642 /** Compute δλ/δv.
643 * @param X satellite position component along f, equinoctial reference frame 1st vector
644 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
645 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
646 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
647 * @return δλ/δv
648 */
649 private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
650 final Vector3D pos = new Vector3D(X, f, Y, g);
651 final Vector3D v2 = new Vector3D(k, getHoV(X, Y, Xdot, Ydot), -h, getKoV(X, Y, Xdot, Ydot));
652 return new Vector3D(-2. * ooA, pos, ooBpo, v2, (I * q * Y - p * X) * ooA, w);
653 }
654
655 }
656
657 /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
658 * a {@link org.hipparchus.analysis.UnivariateVectorFunction function}
659 * of the orbital elements using the Gaussian quadrature rule to get the acceleration.
660 */
661 private static class GaussQuadrature {
662
663 // CHECKSTYLE: stop NoWhitespaceAfter
664
665 // Points and weights for the available quadrature orders
666
667 /** Points for quadrature of order 12. */
668 private static final double[] P_12 = {
669 -0.98156063424671910000,
670 -0.90411725637047490000,
671 -0.76990267419430470000,
672 -0.58731795428661740000,
673 -0.36783149899818024000,
674 -0.12523340851146890000,
675 0.12523340851146890000,
676 0.36783149899818024000,
677 0.58731795428661740000,
678 0.76990267419430470000,
679 0.90411725637047490000,
680 0.98156063424671910000
681 };
682
683 /** Weights for quadrature of order 12. */
684 private static final double[] W_12 = {
685 0.04717533638651220000,
686 0.10693932599531830000,
687 0.16007832854334633000,
688 0.20316742672306584000,
689 0.23349253653835478000,
690 0.24914704581340286000,
691 0.24914704581340286000,
692 0.23349253653835478000,
693 0.20316742672306584000,
694 0.16007832854334633000,
695 0.10693932599531830000,
696 0.04717533638651220000
697 };
698
699 /** Points for quadrature of order 16. */
700 private static final double[] P_16 = {
701 -0.98940093499164990000,
702 -0.94457502307323260000,
703 -0.86563120238783160000,
704 -0.75540440835500310000,
705 -0.61787624440264380000,
706 -0.45801677765722737000,
707 -0.28160355077925890000,
708 -0.09501250983763745000,
709 0.09501250983763745000,
710 0.28160355077925890000,
711 0.45801677765722737000,
712 0.61787624440264380000,
713 0.75540440835500310000,
714 0.86563120238783160000,
715 0.94457502307323260000,
716 0.98940093499164990000
717 };
718
719 /** Weights for quadrature of order 16. */
720 private static final double[] W_16 = {
721 0.02715245941175405800,
722 0.06225352393864777000,
723 0.09515851168249283000,
724 0.12462897125553388000,
725 0.14959598881657685000,
726 0.16915651939500256000,
727 0.18260341504492360000,
728 0.18945061045506847000,
729 0.18945061045506847000,
730 0.18260341504492360000,
731 0.16915651939500256000,
732 0.14959598881657685000,
733 0.12462897125553388000,
734 0.09515851168249283000,
735 0.06225352393864777000,
736 0.02715245941175405800
737 };
738
739 /** Points for quadrature of order 20. */
740 private static final double[] P_20 = {
741 -0.99312859918509490000,
742 -0.96397192727791390000,
743 -0.91223442825132600000,
744 -0.83911697182221890000,
745 -0.74633190646015080000,
746 -0.63605368072651510000,
747 -0.51086700195082700000,
748 -0.37370608871541955000,
749 -0.22778585114164507000,
750 -0.07652652113349734000,
751 0.07652652113349734000,
752 0.22778585114164507000,
753 0.37370608871541955000,
754 0.51086700195082700000,
755 0.63605368072651510000,
756 0.74633190646015080000,
757 0.83911697182221890000,
758 0.91223442825132600000,
759 0.96397192727791390000,
760 0.99312859918509490000
761 };
762
763 /** Weights for quadrature of order 20. */
764 private static final double[] W_20 = {
765 0.01761400713915226400,
766 0.04060142980038684000,
767 0.06267204833410904000,
768 0.08327674157670477000,
769 0.10193011981724048000,
770 0.11819453196151844000,
771 0.13168863844917678000,
772 0.14209610931838212000,
773 0.14917298647260380000,
774 0.15275338713072600000,
775 0.15275338713072600000,
776 0.14917298647260380000,
777 0.14209610931838212000,
778 0.13168863844917678000,
779 0.11819453196151844000,
780 0.10193011981724048000,
781 0.08327674157670477000,
782 0.06267204833410904000,
783 0.04060142980038684000,
784 0.01761400713915226400
785 };
786
787 /** Points for quadrature of order 24. */
788 private static final double[] P_24 = {
789 -0.99518721999702130000,
790 -0.97472855597130950000,
791 -0.93827455200273270000,
792 -0.88641552700440100000,
793 -0.82000198597390300000,
794 -0.74012419157855440000,
795 -0.64809365193697550000,
796 -0.54542147138883950000,
797 -0.43379350762604520000,
798 -0.31504267969616340000,
799 -0.19111886747361634000,
800 -0.06405689286260563000,
801 0.06405689286260563000,
802 0.19111886747361634000,
803 0.31504267969616340000,
804 0.43379350762604520000,
805 0.54542147138883950000,
806 0.64809365193697550000,
807 0.74012419157855440000,
808 0.82000198597390300000,
809 0.88641552700440100000,
810 0.93827455200273270000,
811 0.97472855597130950000,
812 0.99518721999702130000
813 };
814
815 /** Weights for quadrature of order 24. */
816 private static final double[] W_24 = {
817 0.01234122979998733500,
818 0.02853138862893380600,
819 0.04427743881741981000,
820 0.05929858491543691500,
821 0.07334648141108027000,
822 0.08619016153195320000,
823 0.09761865210411391000,
824 0.10744427011596558000,
825 0.11550566805372553000,
826 0.12167047292780335000,
827 0.12583745634682825000,
828 0.12793819534675221000,
829 0.12793819534675221000,
830 0.12583745634682825000,
831 0.12167047292780335000,
832 0.11550566805372553000,
833 0.10744427011596558000,
834 0.09761865210411391000,
835 0.08619016153195320000,
836 0.07334648141108027000,
837 0.05929858491543691500,
838 0.04427743881741981000,
839 0.02853138862893380600,
840 0.01234122979998733500
841 };
842
843 /** Points for quadrature of order 32. */
844 private static final double[] P_32 = {
845 -0.99726386184948160000,
846 -0.98561151154526840000,
847 -0.96476225558750640000,
848 -0.93490607593773970000,
849 -0.89632115576605220000,
850 -0.84936761373256990000,
851 -0.79448379596794250000,
852 -0.73218211874028970000,
853 -0.66304426693021520000,
854 -0.58771575724076230000,
855 -0.50689990893222950000,
856 -0.42135127613063540000,
857 -0.33186860228212767000,
858 -0.23928736225213710000,
859 -0.14447196158279646000,
860 -0.04830766568773831000,
861 0.04830766568773831000,
862 0.14447196158279646000,
863 0.23928736225213710000,
864 0.33186860228212767000,
865 0.42135127613063540000,
866 0.50689990893222950000,
867 0.58771575724076230000,
868 0.66304426693021520000,
869 0.73218211874028970000,
870 0.79448379596794250000,
871 0.84936761373256990000,
872 0.89632115576605220000,
873 0.93490607593773970000,
874 0.96476225558750640000,
875 0.98561151154526840000,
876 0.99726386184948160000
877 };
878
879 /** Weights for quadrature of order 32. */
880 private static final double[] W_32 = {
881 0.00701861000947013600,
882 0.01627439473090571200,
883 0.02539206530926214200,
884 0.03427386291302141000,
885 0.04283589802222658600,
886 0.05099805926237621600,
887 0.05868409347853559000,
888 0.06582222277636193000,
889 0.07234579410884862000,
890 0.07819389578707042000,
891 0.08331192422694673000,
892 0.08765209300440380000,
893 0.09117387869576390000,
894 0.09384439908080441000,
895 0.09563872007927487000,
896 0.09654008851472784000,
897 0.09654008851472784000,
898 0.09563872007927487000,
899 0.09384439908080441000,
900 0.09117387869576390000,
901 0.08765209300440380000,
902 0.08331192422694673000,
903 0.07819389578707042000,
904 0.07234579410884862000,
905 0.06582222277636193000,
906 0.05868409347853559000,
907 0.05099805926237621600,
908 0.04283589802222658600,
909 0.03427386291302141000,
910 0.02539206530926214200,
911 0.01627439473090571200,
912 0.00701861000947013600
913 };
914
915 /** Points for quadrature of order 40. */
916 private static final double[] P_40 = {
917 -0.99823770971055930000,
918 -0.99072623869945710000,
919 -0.97725994998377420000,
920 -0.95791681921379170000,
921 -0.93281280827867660000,
922 -0.90209880696887420000,
923 -0.86595950321225960000,
924 -0.82461223083331170000,
925 -0.77830565142651940000,
926 -0.72731825518992710000,
927 -0.67195668461417960000,
928 -0.61255388966798030000,
929 -0.54946712509512820000,
930 -0.48307580168617870000,
931 -0.41377920437160500000,
932 -0.34199409082575850000,
933 -0.26815218500725370000,
934 -0.19269758070137110000,
935 -0.11608407067525522000,
936 -0.03877241750605081600,
937 0.03877241750605081600,
938 0.11608407067525522000,
939 0.19269758070137110000,
940 0.26815218500725370000,
941 0.34199409082575850000,
942 0.41377920437160500000,
943 0.48307580168617870000,
944 0.54946712509512820000,
945 0.61255388966798030000,
946 0.67195668461417960000,
947 0.72731825518992710000,
948 0.77830565142651940000,
949 0.82461223083331170000,
950 0.86595950321225960000,
951 0.90209880696887420000,
952 0.93281280827867660000,
953 0.95791681921379170000,
954 0.97725994998377420000,
955 0.99072623869945710000,
956 0.99823770971055930000
957 };
958
959 /** Weights for quadrature of order 40. */
960 private static final double[] W_40 = {
961 0.00452127709853309800,
962 0.01049828453115270400,
963 0.01642105838190797300,
964 0.02224584919416689000,
965 0.02793700698002338000,
966 0.03346019528254786500,
967 0.03878216797447199000,
968 0.04387090818567333000,
969 0.04869580763507221000,
970 0.05322784698393679000,
971 0.05743976909939157000,
972 0.06130624249292891000,
973 0.06480401345660108000,
974 0.06791204581523394000,
975 0.07061164739128681000,
976 0.07288658239580408000,
977 0.07472316905796833000,
978 0.07611036190062619000,
979 0.07703981816424793000,
980 0.07750594797842482000,
981 0.07750594797842482000,
982 0.07703981816424793000,
983 0.07611036190062619000,
984 0.07472316905796833000,
985 0.07288658239580408000,
986 0.07061164739128681000,
987 0.06791204581523394000,
988 0.06480401345660108000,
989 0.06130624249292891000,
990 0.05743976909939157000,
991 0.05322784698393679000,
992 0.04869580763507221000,
993 0.04387090818567333000,
994 0.03878216797447199000,
995 0.03346019528254786500,
996 0.02793700698002338000,
997 0.02224584919416689000,
998 0.01642105838190797300,
999 0.01049828453115270400,
1000 0.00452127709853309800
1001 };
1002
1003 /** Points for quadrature of order 48. */
1004 private static final double[] P_48 = {
1005 -0.99877100725242610000,
1006 -0.99353017226635080000,
1007 -0.98412458372282700000,
1008 -0.97059159254624720000,
1009 -0.95298770316043080000,
1010 -0.93138669070655440000,
1011 -0.90587913671556960000,
1012 -0.87657202027424800000,
1013 -0.84358826162439350000,
1014 -0.80706620402944250000,
1015 -0.76715903251574020000,
1016 -0.72403413092381470000,
1017 -0.67787237963266400000,
1018 -0.62886739677651370000,
1019 -0.57722472608397270000,
1020 -0.52316097472223300000,
1021 -0.46690290475095840000,
1022 -0.40868648199071680000,
1023 -0.34875588629216070000,
1024 -0.28736248735545555000,
1025 -0.22476379039468908000,
1026 -0.16122235606889174000,
1027 -0.09700469920946270000,
1028 -0.03238017096286937000,
1029 0.03238017096286937000,
1030 0.09700469920946270000,
1031 0.16122235606889174000,
1032 0.22476379039468908000,
1033 0.28736248735545555000,
1034 0.34875588629216070000,
1035 0.40868648199071680000,
1036 0.46690290475095840000,
1037 0.52316097472223300000,
1038 0.57722472608397270000,
1039 0.62886739677651370000,
1040 0.67787237963266400000,
1041 0.72403413092381470000,
1042 0.76715903251574020000,
1043 0.80706620402944250000,
1044 0.84358826162439350000,
1045 0.87657202027424800000,
1046 0.90587913671556960000,
1047 0.93138669070655440000,
1048 0.95298770316043080000,
1049 0.97059159254624720000,
1050 0.98412458372282700000,
1051 0.99353017226635080000,
1052 0.99877100725242610000
1053 };
1054
1055 /** Weights for quadrature of order 48. */
1056 private static final double[] W_48 = {
1057 0.00315334605230596250,
1058 0.00732755390127620800,
1059 0.01147723457923446900,
1060 0.01557931572294386600,
1061 0.01961616045735556700,
1062 0.02357076083932435600,
1063 0.02742650970835688000,
1064 0.03116722783279807000,
1065 0.03477722256477045000,
1066 0.03824135106583080600,
1067 0.04154508294346483000,
1068 0.04467456085669424000,
1069 0.04761665849249054000,
1070 0.05035903555385448000,
1071 0.05289018948519365000,
1072 0.05519950369998416500,
1073 0.05727729210040315000,
1074 0.05911483969839566000,
1075 0.06070443916589384000,
1076 0.06203942315989268000,
1077 0.06311419228625403000,
1078 0.06392423858464817000,
1079 0.06446616443595010000,
1080 0.06473769681268386000,
1081 0.06473769681268386000,
1082 0.06446616443595010000,
1083 0.06392423858464817000,
1084 0.06311419228625403000,
1085 0.06203942315989268000,
1086 0.06070443916589384000,
1087 0.05911483969839566000,
1088 0.05727729210040315000,
1089 0.05519950369998416500,
1090 0.05289018948519365000,
1091 0.05035903555385448000,
1092 0.04761665849249054000,
1093 0.04467456085669424000,
1094 0.04154508294346483000,
1095 0.03824135106583080600,
1096 0.03477722256477045000,
1097 0.03116722783279807000,
1098 0.02742650970835688000,
1099 0.02357076083932435600,
1100 0.01961616045735556700,
1101 0.01557931572294386600,
1102 0.01147723457923446900,
1103 0.00732755390127620800,
1104 0.00315334605230596250
1105 };
1106 // CHECKSTYLE: resume NoWhitespaceAfter
1107
1108 /** Node points. */
1109 private final double[] nodePoints;
1110
1111 /** Node weights. */
1112 private final double[] nodeWeights;
1113
1114 /** Creates a Gauss integrator of the given order.
1115 *
1116 * @param numberOfPoints Order of the integration rule.
1117 */
1118 GaussQuadrature(final int numberOfPoints) {
1119
1120 switch(numberOfPoints) {
1121 case 12 :
1122 this.nodePoints = P_12.clone();
1123 this.nodeWeights = W_12.clone();
1124 break;
1125 case 16 :
1126 this.nodePoints = P_16.clone();
1127 this.nodeWeights = W_16.clone();
1128 break;
1129 case 20 :
1130 this.nodePoints = P_20.clone();
1131 this.nodeWeights = W_20.clone();
1132 break;
1133 case 24 :
1134 this.nodePoints = P_24.clone();
1135 this.nodeWeights = W_24.clone();
1136 break;
1137 case 32 :
1138 this.nodePoints = P_32.clone();
1139 this.nodeWeights = W_32.clone();
1140 break;
1141 case 40 :
1142 this.nodePoints = P_40.clone();
1143 this.nodeWeights = W_40.clone();
1144 break;
1145 case 48 :
1146 default :
1147 this.nodePoints = P_48.clone();
1148 this.nodeWeights = W_48.clone();
1149 break;
1150 }
1151
1152 }
1153
1154 /** Integrates a given function on the given interval.
1155 *
1156 * @param f Function to integrate.
1157 * @param lowerBound Lower bound of the integration interval.
1158 * @param upperBound Upper bound of the integration interval.
1159 * @return the integral of the weighted function.
1160 */
1161 public double[] integrate(final UnivariateVectorFunction f,
1162 final double lowerBound, final double upperBound) {
1163
1164 final double[] adaptedPoints = nodePoints.clone();
1165 final double[] adaptedWeights = nodeWeights.clone();
1166 transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
1167 return basicIntegrate(f, adaptedPoints, adaptedWeights);
1168 }
1169
1170 /** Performs a change of variable so that the integration
1171 * can be performed on an arbitrary interval {@code [a, b]}.
1172 * <p>
1173 * It is assumed that the natural interval is {@code [-1, 1]}.
1174 * </p>
1175 *
1176 * @param points Points to adapt to the new interval.
1177 * @param weights Weights to adapt to the new interval.
1178 * @param a Lower bound of the integration interval.
1179 * @param b Lower bound of the integration interval.
1180 */
1181 private void transform(final double[] points, final double[] weights,
1182 final double a, final double b) {
1183 // Scaling
1184 final double scale = (b - a) / 2;
1185 final double shift = a + scale;
1186 for (int i = 0; i < points.length; i++) {
1187 points[i] = points[i] * scale + shift;
1188 weights[i] *= scale;
1189 }
1190 }
1191
1192 /** Returns an estimate of the integral of {@code f(x) * w(x)},
1193 * where {@code w} is a weight function that depends on the actual
1194 * flavor of the Gauss integration scheme.
1195 *
1196 * @param f Function to integrate.
1197 * @param points Nodes.
1198 * @param weights Nodes weights.
1199 * @return the integral of the weighted function.
1200 */
1201 private double[] basicIntegrate(final UnivariateVectorFunction f,
1202 final double[] points,
1203 final double[] weights) {
1204 double x = points[0];
1205 double w = weights[0];
1206 double[] v = f.value(x);
1207 final double[] y = new double[v.length];
1208 for (int j = 0; j < v.length; j++) {
1209 y[j] = w * v[j];
1210 }
1211 final double[] t = y.clone();
1212 final double[] c = new double[v.length];
1213 final double[] s = t.clone();
1214 for (int i = 1; i < points.length; i++) {
1215 x = points[i];
1216 w = weights[i];
1217 v = f.value(x);
1218 for (int j = 0; j < v.length; j++) {
1219 y[j] = w * v[j] - c[j];
1220 t[j] = s[j] + y[j];
1221 c[j] = (t[j] - s[j]) - y[j];
1222 s[j] = t[j];
1223 }
1224 }
1225 return s;
1226 }
1227
1228 }
1229
1230 /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
1231 * <p>
1232 * Those coefficients are given in Danielson paper by expression 4.4-(6)
1233 * </p>
1234 * @author Petre Bazavan
1235 * @author Lucian Barbulescu
1236 */
1237 private class FourierCjSjCoefficients {
1238
1239 /** Maximum possible value for j. */
1240 private final int jMax;
1241
1242 /** The C<sub>i</sub><sup>j</sup> coefficients.
1243 * <p>
1244 * the index i corresponds to the following elements: <br/>
1245 * - 0 for a <br>
1246 * - 1 for k <br>
1247 * - 2 for h <br>
1248 * - 3 for q <br>
1249 * - 4 for p <br>
1250 * - 5 for λ <br>
1251 * </p>
1252 */
1253 private final double[][] cCoef;
1254
1255 /** The C<sub>i</sub><sup>j</sup> coefficients.
1256 * <p>
1257 * the index i corresponds to the following elements: <br/>
1258 * - 0 for a <br>
1259 * - 1 for k <br>
1260 * - 2 for h <br>
1261 * - 3 for q <br>
1262 * - 4 for p <br>
1263 * - 5 for λ <br>
1264 * </p>
1265 */
1266 private final double[][] sCoef;
1267
1268 /** Standard constructor.
1269 * @param state the current state
1270 * @param jMax maximum value for j
1271 * @throws OrekitException in case of an error
1272 */
1273 FourierCjSjCoefficients(final SpacecraftState state, final int jMax)
1274 throws OrekitException {
1275 //Initialise the fields
1276 this.jMax = jMax;
1277
1278 //Allocate the arrays
1279 final int rows = jMax + 1;
1280 cCoef = new double[rows][6];
1281 sCoef = new double[rows][6];
1282
1283 //Compute the coefficients
1284 computeCoefficients(state);
1285 }
1286
1287 /**
1288 * Compute the Fourrier coefficients.
1289 * <p>
1290 * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
1291 * as D<sub>i</sub><sup>m</sup> is always 0.
1292 * </p>
1293 * @param state the current state
1294 * @throws OrekitException in case of an error
1295 */
1296 private void computeCoefficients(final SpacecraftState state)
1297 throws OrekitException {
1298 // Computes the limits for the integral
1299 final double[] ll = getLLimits(state);
1300 // Computes integrated mean element rates if Llow < Lhigh
1301 if (ll[0] < ll[1]) {
1302 //Compute 1 / PI
1303 final double ooPI = 1 / FastMath.PI;
1304
1305 // loop through all values of j
1306 for (int j = 0; j <= jMax; j++) {
1307 final double[] curentCoefficients =
1308 integrator.integrate(new IntegrableFunction(state, false, j), ll[0], ll[1]);
1309
1310 //divide by PI and set the values for the coefficients
1311 for (int i = 0; i < 6; i++) {
1312 cCoef[j][i] = ooPI * curentCoefficients[i];
1313 sCoef[j][i] = ooPI * curentCoefficients[i + 6];
1314 }
1315 }
1316 }
1317 }
1318
1319 /** Get the coefficient C<sub>i</sub><sup>j</sup>.
1320 * @param i i index - corresponds to the required variation
1321 * @param j j index
1322 * @return the coefficient C<sub>i</sub><sup>j</sup>
1323 */
1324 public double getCij(final int i, final int j) {
1325 return cCoef[j][i];
1326 }
1327
1328 /** Get the coefficient S<sub>i</sub><sup>j</sup>.
1329 * @param i i index - corresponds to the required variation
1330 * @param j j index
1331 * @return the coefficient S<sub>i</sub><sup>j</sup>
1332 */
1333 public double getSij(final int i, final int j) {
1334 return sCoef[j][i];
1335 }
1336 }
1337
1338 /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
1339 *
1340 * <p>
1341 * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
1342 * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m > 2 are also 0.
1343 * </p>
1344 * @author Petre Bazavan
1345 * @author Lucian Barbulescu
1346 *
1347 */
1348 private static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {
1349
1350 /** Serializable UID. */
1351 private static final long serialVersionUID = 20151118L;
1352
1353 /** Maximum value for j index. */
1354 private final int jMax;
1355
1356 /** Number of points used in the interpolation process. */
1357 private final int interpolationPoints;
1358
1359 /** Prefix for coefficients keys. */
1360 private final String coefficientsKeyPrefix;
1361
1362 /** All coefficients slots. */
1363 private final transient TimeSpanMap<Slot> slots;
1364
1365 /** Constructor.
1366 * @param coefficientsKeyPrefix prefix for coefficients keys
1367 * @param jMax maximum value for j index
1368 * @param interpolationPoints number of points used in the interpolation process
1369 * @param slots all coefficients slots
1370 */
1371 GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
1372 final int jMax, final int interpolationPoints,
1373 final TimeSpanMap<Slot> slots) {
1374 //Initialize fields
1375 this.jMax = jMax;
1376 this.interpolationPoints = interpolationPoints;
1377 this.coefficientsKeyPrefix = coefficientsKeyPrefix;
1378 this.slots = slots;
1379 }
1380
1381 /** Get the slot valid for some date.
1382 * @param meanStates mean states defining the slot
1383 * @return slot valid at the specified date
1384 */
1385 public Slot createSlot(final SpacecraftState ... meanStates) {
1386 final Slot slot = new Slot(jMax, interpolationPoints);
1387 final AbsoluteDate first = meanStates[0].getDate();
1388 final AbsoluteDate last = meanStates[meanStates.length - 1].getDate();
1389 if (first.compareTo(last) <= 0) {
1390 slots.addValidAfter(slot, first);
1391 } else {
1392 slots.addValidBefore(slot, first);
1393 }
1394 return slot;
1395 }
1396
1397 /** Compute the short periodic coefficients.
1398 *
1399 * @param state current state information: date, kinematics, attitude
1400 * @param slot coefficients slot
1401 * @param fourierCjSj Fourier coefficients
1402 * @param uijvij U and V coefficients
1403 * @param n Keplerian mean motion
1404 * @param a semi major axis
1405 * @throws OrekitException if an error occurs
1406 */
1407 private void computeCoefficients(final SpacecraftState state, final Slot slot,
1408 final FourierCjSjCoefficients fourierCjSj,
1409 final UijVijCoefficients uijvij,
1410 final double n, final double a)
1411 throws OrekitException {
1412
1413 // get the current date
1414 final AbsoluteDate date = state.getDate();
1415
1416 // compute the k₂⁰ coefficient
1417 final double k20 = computeK20(jMax, uijvij.currentRhoSigmaj);
1418
1419 // 1. / n
1420 final double oon = 1. / n;
1421 // 3. / (2 * a * n)
1422 final double to2an = 1.5 * oon / a;
1423 // 3. / (4 * a * n)
1424 final double to4an = to2an / 2;
1425
1426 // Compute the coefficients for each element
1427 final int size = jMax + 1;
1428 final double[] di1 = new double[6];
1429 final double[] di2 = new double[6];
1430 final double[][] currentCij = new double[size][6];
1431 final double[][] currentSij = new double[size][6];
1432 for (int i = 0; i < 6; i++) {
1433
1434 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
1435 di1[i] = -oon * fourierCjSj.getCij(i, 0);
1436 if (i == 5) {
1437 di1[i] += to2an * uijvij.getU1(0, 0);
1438 }
1439 di2[i] = 0.;
1440 if (i == 5) {
1441 di2[i] += -to4an * fourierCjSj.getCij(0, 0);
1442 }
1443
1444 //the C<sub>i</sub>⁰ is computed based on all others
1445 currentCij[0][i] = -di2[i] * k20;
1446
1447 for (int j = 1; j <= jMax; j++) {
1448 // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
1449 currentCij[j][i] = oon * uijvij.getU1(j, i);
1450 if (i == 5) {
1451 currentCij[j][i] += -to2an * uijvij.getU2(j);
1452 }
1453 currentSij[j][i] = oon * uijvij.getV1(j, i);
1454 if (i == 5) {
1455 currentSij[j][i] += -to2an * uijvij.getV2(j);
1456 }
1457
1458 // add the computed coefficients to C<sub>i</sub>⁰
1459 currentCij[0][i] += -(currentCij[j][i] * uijvij.currentRhoSigmaj[0][j] + currentSij[j][i] * uijvij.currentRhoSigmaj[1][j]);
1460 }
1461
1462 }
1463
1464 // add the values to the interpolators
1465 slot.cij[0].addGridPoint(date, currentCij[0]);
1466 slot.dij[1].addGridPoint(date, di1);
1467 slot.dij[2].addGridPoint(date, di2);
1468 for (int j = 1; j <= jMax; j++) {
1469 slot.cij[j].addGridPoint(date, currentCij[j]);
1470 slot.sij[j].addGridPoint(date, currentSij[j]);
1471 }
1472
1473 }
1474
1475 /** Compute the coefficient k₂⁰ by using the equation
1476 * 2.5.3-(9a) from Danielson.
1477 * <p>
1478 * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
1479 * k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
1480 * </p>
1481 * @param kMax max value fot k index
1482 * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
1483 * @return the coefficient k₂⁰
1484 */
1485 private double computeK20(final int kMax, final double[][] currentRhoSigmaj) {
1486 double k20 = 0.;
1487
1488 for (int kIndex = 1; kIndex <= kMax; kIndex++) {
1489 // After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
1490 //k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
1491 double currentTerm = currentRhoSigmaj[1][kIndex] * currentRhoSigmaj[1][kIndex] +
1492 currentRhoSigmaj[0][kIndex] * currentRhoSigmaj[0][kIndex];
1493
1494 //multiply by 2 / k²
1495 currentTerm *= 2. / (kIndex * kIndex);
1496
1497 // add the term to the result
1498 k20 += currentTerm;
1499 }
1500
1501 return k20;
1502 }
1503
1504 /** {@inheritDoc} */
1505 @Override
1506 public double[] value(final Orbit meanOrbit) {
1507
1508 // select the coefficients slot
1509 final Slot slot = slots.get(meanOrbit.getDate());
1510
1511 // Get the True longitude L
1512 final double L = meanOrbit.getLv();
1513
1514 // Compute the center (l - λ)
1515 final double center = L - meanOrbit.getLM();
1516 // Compute (l - λ)²
1517 final double center2 = center * center;
1518
1519 // Initialize short periodic variations
1520 final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
1521 final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
1522 final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
1523 for (int i = 0; i < 6; i++) {
1524 shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
1525 }
1526
1527 for (int j = 1; j <= JMAX; j++) {
1528 final double[] c = slot.cij[j].value(meanOrbit.getDate());
1529 final double[] s = slot.sij[j].value(meanOrbit.getDate());
1530 final double cos = FastMath.cos(j * L);
1531 final double sin = FastMath.sin(j * L);
1532 for (int i = 0; i < 6; i++) {
1533 // add corresponding term to the short periodic variation
1534 shortPeriodicVariation[i] += c[i] * cos;
1535 shortPeriodicVariation[i] += s[i] * sin;
1536 }
1537 }
1538
1539 return shortPeriodicVariation;
1540
1541 }
1542
1543 /** {@inheritDoc} */
1544 public String getCoefficientsKeyPrefix() {
1545 return coefficientsKeyPrefix;
1546 }
1547
1548 /** {@inheritDoc}
1549 * <p>
1550 * For Gaussian forces,there are JMAX cj coefficients,
1551 * JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
1552 * this sums up to 27 coefficients. The j index is the integer
1553 * multiplier for the true longitude argument in the cj and sj
1554 * coefficients and to the degree in the polynomial dj coefficients.
1555 * </p>
1556 */
1557 @Override
1558 public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
1559 throws OrekitException {
1560
1561 // select the coefficients slot
1562 final Slot slot = slots.get(date);
1563
1564 final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * JMAX + 3);
1565 storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
1566 storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
1567 storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
1568 for (int j = 1; j <= JMAX; j++) {
1569 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
1570 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
1571 }
1572
1573 return coefficients;
1574
1575 }
1576
1577 /** Put a coefficient in a map if selected.
1578 * @param map map to populate
1579 * @param selected set of coefficients that should be put in the map
1580 * (empty set means all coefficients are selected)
1581 * @param value coefficient value
1582 * @param id coefficient identifier
1583 * @param indices list of coefficient indices
1584 */
1585 private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
1586 final double[] value, final String id, final int ... indices) {
1587 final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
1588 keyBuilder.append(id);
1589 for (int index : indices) {
1590 keyBuilder.append('[').append(index).append(']');
1591 }
1592 final String key = keyBuilder.toString();
1593 if (selected.isEmpty() || selected.contains(key)) {
1594 map.put(key, value);
1595 }
1596 }
1597
1598 /** Replace the instance with a data transfer object for serialization.
1599 * @return data transfer object that will be serialized
1600 * @exception NotSerializableException if an additional state provider is not serializable
1601 */
1602 private Object writeReplace() throws NotSerializableException {
1603
1604 // slots transitions
1605 final SortedSet<TimeSpanMap.Transition<Slot>> transitions = slots.getTransitions();
1606 final AbsoluteDate[] transitionDates = new AbsoluteDate[transitions.size()];
1607 final Slot[] allSlots = new Slot[transitions.size() + 1];
1608 int i = 0;
1609 for (final TimeSpanMap.Transition<Slot> transition : transitions) {
1610 if (i == 0) {
1611 // slot before the first transition
1612 allSlots[i] = transition.getBefore();
1613 }
1614 if (i < transitionDates.length) {
1615 transitionDates[i] = transition.getDate();
1616 allSlots[++i] = transition.getAfter();
1617 }
1618 }
1619
1620 return new DataTransferObject(jMax, interpolationPoints, coefficientsKeyPrefix,
1621 transitionDates, allSlots);
1622
1623 }
1624
1625
1626 /** Internal class used only for serialization. */
1627 private static class DataTransferObject implements Serializable {
1628
1629 /** Serializable UID. */
1630 private static final long serialVersionUID = 20160319L;
1631
1632 /** Maximum value for j index. */
1633 private final int jMax;
1634
1635 /** Number of points used in the interpolation process. */
1636 private final int interpolationPoints;
1637
1638 /** Prefix for coefficients keys. */
1639 private final String coefficientsKeyPrefix;
1640
1641 /** Transitions dates. */
1642 private final AbsoluteDate[] transitionDates;
1643
1644 /** All slots. */
1645 private final Slot[] allSlots;
1646
1647 /** Simple constructor.
1648 * @param jMax maximum value for j index
1649 * @param interpolationPoints number of points used in the interpolation process
1650 * @param coefficientsKeyPrefix prefix for coefficients keys
1651 * @param transitionDates transitions dates
1652 * @param allSlots all slots
1653 */
1654 DataTransferObject(final int jMax, final int interpolationPoints,
1655 final String coefficientsKeyPrefix,
1656 final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
1657 this.jMax = jMax;
1658 this.interpolationPoints = interpolationPoints;
1659 this.coefficientsKeyPrefix = coefficientsKeyPrefix;
1660 this.transitionDates = transitionDates;
1661 this.allSlots = allSlots;
1662 }
1663
1664 /** Replace the deserialized data transfer object with a {@link GaussianShortPeriodicCoefficients}.
1665 * @return replacement {@link GaussianShortPeriodicCoefficients}
1666 */
1667 private Object readResolve() {
1668
1669 final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
1670 for (int i = 0; i < transitionDates.length; ++i) {
1671 slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
1672 }
1673
1674 return new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix, jMax, interpolationPoints, slots);
1675
1676 }
1677
1678 }
1679
1680 }
1681
1682 /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
1683 * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
1684 * <p>
1685 * The index i takes only the values 1 and 2<br>
1686 * For U only the index 0 for j is used.
1687 * </p>
1688 *
1689 * @author Petre Bazavan
1690 * @author Lucian Barbulescu
1691 */
1692 private static class UijVijCoefficients {
1693
1694 /** The U₁<sup>j</sup> coefficients.
1695 * <p>
1696 * The first index identifies the Fourier coefficients used<br>
1697 * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
1698 * The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
1699 * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
1700 * to compute the coefficients U₂<sup>j</sup>
1701 * </p>
1702 */
1703 private final double[][] u1ij;
1704
1705 /** The V₁<sup>j</sup> coefficients.
1706 * <p>
1707 * The first index identifies the Fourier coefficients used<br>
1708 * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
1709 * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
1710 * to compute the coefficients V₂<sup>j</sup>
1711 * </p>
1712 */
1713 private final double[][] v1ij;
1714
1715 /** The U₂<sup>j</sup> coefficients.
1716 * <p>
1717 * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
1718 * </p>
1719 */
1720 private final double[] u2ij;
1721
1722 /** The V₂<sup>j</sup> coefficients.
1723 * <p>
1724 * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
1725 * </p>
1726 */
1727 private final double[] v2ij;
1728
1729 /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
1730 private final double[][] currentRhoSigmaj;
1731
1732 /** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
1733 private final FourierCjSjCoefficients fourierCjSj;
1734
1735 /** The maximum value for j index. */
1736 private final int jMax;
1737
1738 /** Constructor.
1739 * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
1740 * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
1741 * @param jMax maximum value for j index
1742 */
1743 UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
1744 this.currentRhoSigmaj = currentRhoSigmaj;
1745 this.fourierCjSj = fourierCjSj;
1746 this.jMax = jMax;
1747
1748 // initialize the internal arrays.
1749 this.u1ij = new double[6][2 * jMax + 1];
1750 this.v1ij = new double[6][2 * jMax + 1];
1751 this.u2ij = new double[jMax + 1];
1752 this.v2ij = new double[jMax + 1];
1753
1754 //compute the coefficients
1755 computeU1V1Coefficients();
1756 computeU2V2Coefficients();
1757 }
1758
1759 /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
1760 private void computeU1V1Coefficients() {
1761 // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
1762 // for j >= 1
1763 // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
1764 u1ij[0][0] = 0;
1765 for (int j = 1; j <= jMax; j++) {
1766 // compute 1 / j
1767 final double ooj = 1. / j;
1768
1769 for (int i = 0; i < 6; i++) {
1770 //j is aready between 1 and J
1771 u1ij[i][j] = fourierCjSj.getSij(i, j);
1772 v1ij[i][j] = fourierCjSj.getCij(i, j);
1773
1774 // 1 - δ<sub>1j</sub> is 1 for all j > 1
1775 if (j > 1) {
1776 // k starts with 1 because j-J is less than or equal to 0
1777 for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
1778 // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
1779 // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
1780 u1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
1781 fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];
1782
1783 // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
1784 // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
1785 v1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
1786 fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
1787 }
1788 }
1789
1790 // since j must be between 1 and J-1 and is already between 1 and J
1791 // the following sum is skiped only for j = jMax
1792 if (j != jMax) {
1793 for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
1794 // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
1795 // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
1796 u1ij[i][j] += -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
1797 fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];
1798
1799 // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
1800 // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
1801 v1ij[i][j] += fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
1802 fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
1803 }
1804 }
1805
1806 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
1807 // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
1808 // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
1809 u1ij[i][j] += -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
1810 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];
1811
1812 // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
1813 // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
1814 v1ij[i][j] += fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
1815 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
1816 }
1817
1818 // divide by 1 / j
1819 u1ij[i][j] *= -ooj;
1820 v1ij[i][j] *= ooj;
1821
1822 // if index = 1 (i == 0) add the computed terms to U₁⁰
1823 if (i == 0) {
1824 //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
1825 u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
1826 }
1827 }
1828 }
1829
1830 // Terms with j > jMax are required only when computing the coefficients
1831 // U₂<sup>j</sup> and V₂<sup>j</sup>
1832 // and those coefficients are only required for Fourier index = 1 (i == 0).
1833 for (int j = jMax + 1; j <= 2 * jMax; j++) {
1834 // compute 1 / j
1835 final double ooj = 1. / j;
1836 //the value of i is 0
1837 u1ij[0][j] = 0.;
1838 v1ij[0][j] = 0.;
1839
1840 //k starts from j-J as it is always greater than or equal to 1
1841 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
1842 // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
1843 // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
1844 u1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
1845 fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];
1846
1847 // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
1848 // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
1849 v1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
1850 fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
1851 }
1852 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
1853 // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
1854 // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
1855 u1ij[0][j] += -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
1856 fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];
1857
1858 // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
1859 // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
1860 v1ij[0][j] += fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
1861 fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
1862 }
1863
1864 // divide by 1 / j
1865 u1ij[0][j] *= -ooj;
1866 v1ij[0][j] *= ooj;
1867 }
1868 }
1869
1870 /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
1871 * <p>
1872 * Only the coefficients for Fourier index = 1 (i == 0) are required.
1873 * </p>
1874 */
1875 private void computeU2V2Coefficients() {
1876 for (int j = 1; j <= jMax; j++) {
1877 // compute 1 / j
1878 final double ooj = 1. / j;
1879
1880 // only the values for i == 0 are computed
1881 u2ij[j] = v1ij[0][j];
1882 v2ij[j] = u1ij[0][j];
1883
1884 // 1 - δ<sub>1j</sub> is 1 for all j > 1
1885 if (j > 1) {
1886 for (int l = 1; l <= j - 1; l++) {
1887 // U₁<sup>j-l</sup> * σ<sub>l</sub> +
1888 // V₁<sup>j-l</sup> * ρ<sub>l</sub>
1889 u2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[1][l] +
1890 v1ij[0][j - l] * currentRhoSigmaj[0][l];
1891
1892 // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
1893 // V₁<sup>j-l</sup> * σ<sub>l</sub>
1894 v2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[0][l] -
1895 v1ij[0][j - l] * currentRhoSigmaj[1][l];
1896 }
1897 }
1898
1899 for (int l = 1; l <= jMax; l++) {
1900 // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
1901 // U₁<sup>l</sup> * σ<sub>j+l</sub> +
1902 // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
1903 // V₁<sup>l</sup> * ρ<sub>j+l</sub>
1904 u2ij[j] += -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
1905 u1ij[0][l] * currentRhoSigmaj[1][j + l] +
1906 v1ij[0][j + l] * currentRhoSigmaj[0][l] -
1907 v1ij[0][l] * currentRhoSigmaj[0][j + l];
1908
1909 // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
1910 // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
1911 // V₁<sup>j+l</sup> * σ<sub>l</sub> +
1912 // V₁<sup>l</sup> * σ<sub>j+l</sub>
1913 u2ij[j] += u1ij[0][j + l] * currentRhoSigmaj[0][l] +
1914 u1ij[0][l] * currentRhoSigmaj[0][j + l] +
1915 v1ij[0][j + l] * currentRhoSigmaj[1][l] +
1916 v1ij[0][l] * currentRhoSigmaj[1][j + l];
1917 }
1918
1919 // divide by 1 / j
1920 u2ij[j] *= -ooj;
1921 v2ij[j] *= ooj;
1922 }
1923 }
1924
1925 /** Get the coefficient U₁<sup>j</sup> for Fourier index i.
1926 *
1927 * @param j j index
1928 * @param i Fourier index (starts at 0)
1929 * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
1930 */
1931 public double getU1(final int j, final int i) {
1932 return u1ij[i][j];
1933 }
1934
1935 /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
1936 *
1937 * @param j j index
1938 * @param i Fourier index (starts at 0)
1939 * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
1940 */
1941 public double getV1(final int j, final int i) {
1942 return v1ij[i][j];
1943 }
1944
1945 /** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
1946 *
1947 * @param j j index
1948 * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
1949 */
1950 public double getU2(final int j) {
1951 return u2ij[j];
1952 }
1953
1954 /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
1955 *
1956 * @param j j index
1957 * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
1958 */
1959 public double getV2(final int j) {
1960 return v2ij[j];
1961 }
1962 }
1963
1964 /** Coefficients valid for one time slot. */
1965 private static class Slot implements Serializable {
1966
1967 /** Serializable UID. */
1968 private static final long serialVersionUID = 20160319L;
1969
1970 /**The coefficients D<sub>i</sub><sup>j</sup>.
1971 * <p>
1972 * Only for j = 1 and j = 2 the coefficients are not 0. <br>
1973 * i corresponds to the equinoctial element, as follows:
1974 * - i=0 for a <br/>
1975 * - i=1 for k <br/>
1976 * - i=2 for h <br/>
1977 * - i=3 for q <br/>
1978 * - i=4 for p <br/>
1979 * - i=5 for λ <br/>
1980 * </p>
1981 */
1982 private final ShortPeriodicsInterpolatedCoefficient[] dij;
1983
1984 /** The coefficients C<sub>i</sub><sup>j</sup>.
1985 * <p>
1986 * The index order is cij[j][i] <br/>
1987 * i corresponds to the equinoctial element, as follows: <br/>
1988 * - i=0 for a <br/>
1989 * - i=1 for k <br/>
1990 * - i=2 for h <br/>
1991 * - i=3 for q <br/>
1992 * - i=4 for p <br/>
1993 * - i=5 for λ <br/>
1994 * </p>
1995 */
1996 private final ShortPeriodicsInterpolatedCoefficient[] cij;
1997
1998 /** The coefficients S<sub>i</sub><sup>j</sup>.
1999 * <p>
2000 * The index order is sij[j][i] <br/>
2001 * i corresponds to the equinoctial element, as follows: <br/>
2002 * - i=0 for a <br/>
2003 * - i=1 for k <br/>
2004 * - i=2 for h <br/>
2005 * - i=3 for q <br/>
2006 * - i=4 for p <br/>
2007 * - i=5 for λ <br/>
2008 * </p>
2009 */
2010 private final ShortPeriodicsInterpolatedCoefficient[] sij;
2011
2012 /** Simple constructor.
2013 * @param jMax maximum value for j index
2014 * @param interpolationPoints number of points used in the interpolation process
2015 */
2016 Slot(final int jMax, final int interpolationPoints) {
2017
2018 dij = new ShortPeriodicsInterpolatedCoefficient[3];
2019 cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
2020 sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
2021
2022 // Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
2023 for (int j = 0; j <= jMax; j++) {
2024 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
2025 if (j > 0) {
2026 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
2027 }
2028 // Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
2029 if (j == 1 || j == 2) {
2030 dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
2031 }
2032 }
2033
2034 }
2035
2036 }
2037
2038 }